import sys
sys.path.insert(0,'..')
import Jason2_py as Jason2
import numpy as np
import pandas as pd
# Plotting Functions
import holoviews as hv
import cartopy.crs as ccrs
import plot_funcs as pf
from IPython.core.display import HTML
The Python package written for this project, Jason2_py, can be found at https://github.com/theo-yang/jason2. The package contains the following modules:
Jason2_py.download.py: tools for downloading and reading Jason-2 ssha data directly from the online repository, as well as functions for reading netcdf files from Copernicus datasets. ECCOv4 datasets can be read using ecco_v4_py
Jason2_py.preprocess.py: tools for removing ssha datasets based on NaN and distance thresholds.
find_2D_power_spectrum.py: function for calculating the Hann-windowed two dimensional Fourier transform of ssha data.
Jason2_py.model.py: tools for modeling the two-dimensional Fourier transform of ssha, as discussed here
See full dosctrings using ? operator
The following analysis uses datasets from the following sources:
Now let's download the data for this demonstration (you can also just download from saved data files):
Jason-2
# get existing cycles:
f = open("cycles.txt", "r")
cycles = [line.strip()[1:4] for line in f]
# define reference cycle
ref_cycle = '100'
# target start point and track length
starts = [[14.5,-167],
[13.5,-164.5],
[12.5,-162],
[10.5,-160]]
track_length = 1000
direction = 'N'
# download and read ssha data for ascending tracks South of Hawaiian Ridge
Hawaii_SA_raw = Jason2.nested_dict()
for track in range(4):
# define start and reference file:
start = starts[track]
reference_file = Jason2.download.find_track(start,ref_cycle,direction = direction)
track_name = 'SA' + str(track + 1)
# download, read data from each cycle
for cycle in cycles:
cycle_data = Jason2.download.read_track(start,track_length,ref_cycle,cycle,reference_file)
if cycle_data == None:
continue
Hawaii_SA_raw[track_name][cycle] = cycle_data
Copernicus
(NOTE: change directory names to where nc files are stored)
# set directories to retrieve nc files, initialize data structure
dir_names = ['D:/copernicus_current_data/{year}/{year}'.format(year=year) for year in np.arange(2008,2017)]
current_data = dict()
Hawaii_SA = Jason2.nested_dict()
for t_num,t_name in enumerate(Hawaii_SA_raw.keys()):
# preprocess data
Hawaii_SA[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SA_raw[t_name],starts[t_num],distance_tol=15,masked_tol=5,supress_message=True)
# find latitude and longitude range
lats = np.array([[Hawaii_SA[t_name][cycle]['lat'][0],Hawaii_SA[t_name][cycle]['lat'][-1]] for cycle in Hawaii_SA[t_name].keys()])
lons = np.array([[Hawaii_SA[t_name][cycle]['lon'][0],Hawaii_SA[t_name][cycle]['lon'][-1]] for cycle in Hawaii_SA[t_name].keys()])
lat_range = [min(lats[:,0]), max(lats[:,1])]
lon_range = [min(lons[:,0]), max(lons[:,1])]
# download datasets from .nc directory (change directory name as necessary)
time,lat,lon,u,v = Jason2.download.read_current_data(lat_range,lon_range,dir_names)
current_data[t_name] = dict(zip(['time','lat','lon','u','v'],[time,lat,lon,u,v]))
ECCO version 4
(1) Use the bespoke Python package found here to read in $\frac{d \rho}{dr}$ datasets.
(2) Solve for buoyancy frequency using Jason2_py.find_N_from_DRHODR()
(3) Use Dedalus to solve Sturm-Liouville BVP for horizontal velocities of internal tides (see derivation). This can be done using script dedalus\solve_vertical_modes.py
For simplicity, we can just load the results from mode_plot.py from saved files:
# Jason-2
Hawaii_SA_raw = np.load('./test_data/Hawaii_SA_raw.npy',allow_pickle='TRUE').item()
# Copernicus
current_data = np.load('./test_data/current_dataSA.npy',allow_pickle='TRUE').item()
# ECCOv4
models = Jason2.nested_dict()
tracks = ['SA' + str(i+1) for i in range(4)]
for track in tracks:
modes = 1 / np.genfromtxt('./test_data/eigenvectors/modes%s.csv'%track,delimiter = ',')[:,1]
models[track]['M2 model']['wavenumber'] = modes[0]
models[track]['S2 model']['wavenumber'] = modes[1]
models[track]['eigenvectors'] = np.genfromtxt('./test_data\eigenvectors\%s.csv'%track,delimiter = ',')
Plot 2D power spectrum
Now that the data is downloaded, let's look compute the 2D power spectrum using Jason2_py.find_2D_power_spectrum():
# initialize plotting
Hawaii_SA = Jason2.nested_dict()
# starting coordinates
starts = [[14.5,-167],
[13.5,-164.5],
[12.5,-162],
[10.5,-160]]
for t_num, t_name in enumerate(Hawaii_SA_raw.keys()):
# preprocess dataset
Hawaii_SA[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SA_raw[t_name],starts[t_num],distance_tol=15,masked_tol=5,supress_message=True)
# find power spectrum
fft_2D, dx, dt, nu_vector, f_vector = Jason2.find_2D_power_spectrum(Hawaii_SA[t_name])
# find average lat and lon range
lats = np.array([[Hawaii_SA[t_name][cycle]['lat'][0],Hawaii_SA[t_name][cycle]['lat'][-1]] for cycle in Hawaii_SA[t_name].keys()])
lons = np.array([[Hawaii_SA[t_name][cycle]['lon'][0],Hawaii_SA[t_name][cycle]['lon'][-1]] for cycle in Hawaii_SA[t_name].keys()])
lat_range = [np.mean(lats[:,0]), np.mean(lats[:,1])]
lon_range = [np.mean(lons[:,0]), np.mean(lons[:,1])]
# save info to models
models[t_name]['data'].update({'fft 2D' : fft_2D,
'dx': dx,
'dt': dt,
'wavenumbers' : nu_vector,
'frequencies' : f_vector,
'lat range' : lat_range,
'lon range' : lon_range})
# draw plots
track_dict = {'Track ' + str(i + 1) : pf.track_plt(models[track]['data'], 'Track ' + str(i + 1)) for i, track in enumerate(tracks)}
power_dict = {'Track ' + str(i + 1) : pf.fft_2D_plt(models[track]['data'], '') for i, track in enumerate(tracks)}
track_plot = hv.HoloMap(track_dict, kdims='Track').opts(infer_projection=True)
power_plot = hv.HoloMap(power_dict, kdims='Track')
(track_plot + power_plot).opts(title='Figure 1: 2D Power Spectra')
Anisotropic Wave
We begin model tides as sinusoids given by the function:
$$h(x,t) = A sin(2 \pi (k x + \omega t))$$where $A$ is wave amplitude, $k$ is the tidal wavenumber and $\omega$ is the alias frequency. In our analysis, we minimize the effects of discontiuities by multiplying by a Hann window in both the space (x) and time (t). This window function is:
$$w(x,t) = \frac{1}{4}(1+cos(\frac{2 \pi x}{L})(1+cos(\frac{2 \pi t}{T}))rect(\frac{x}{L})rect(\frac{t}{T})$$where $L$ is the length of the track and $T$ is the length of the time window. Finally, to account for the discretization of the signal, we multiply by a Dirac comb given by:
$$d(x,t) = \sum^\infty_{n=-\infty}\delta(t - n \Delta t)\sum^\infty_{n=-\infty}\delta(x - n \Delta x)$$where $\Delta t$ and $\Delta x$ are the time and space intervals between sucessive measurements. In total, our signal function is $s(x,t) = h(x,t) \ w(x,t) \ d(x,t)$. We would like to generate an analytical function for the two dimensional Fourier transform, $\hat{s}(\nu,f) = \mathcal{F}[s(x,t)]$. Taking the fourier transforms of $h(x,t)$, $w(x,t)$, and $d(x,t)$, we have:
$$ \begin{align} \hat{h}(\nu,f) &= \frac{1}{2}A i \Big(\delta (\nu + k) \delta (f + \omega) - \delta (\nu - k) \delta (f - \omega) \Big) \\ \hat{w}(\nu,f) &= \frac{1}{4} \Big(\mathcal{F}[1 + cos(\frac{2 \pi x}{L})] * \mathcal{F}[rect(\frac{1}{L})]\Big)\Big(\mathcal{F}[1 + cos(\frac{2 \pi t}{T})] * \mathcal{F}[rect(\frac{t}{T})]\Big) \\ &= \frac{1}{4} \Big[ \Big(\delta(\nu) + \frac{1}{2} \delta (\nu - 1/L) + \frac{1}{2} \delta (\nu + 1/L)\Big) * L sinc(L \nu) \Big] \Big[\Big(\delta(f) + \frac{1}{2} \delta (f - 1/T) + \frac{1}{2} \delta (f + 1/T)\Big) * T sinc(T f) \Big] \\ &= \frac{1}{4}\Big[ L sinc(L \nu) + \frac{1}{2} L sinc(L \nu - 1) + \frac{1}{2} L sinc(L \nu + 1)\Big]\Big[ T sinc(T f) + \frac{1}{2} T sinc(T f - 1) + \frac{1}{2} T sinc(T f + 1)\Big] \\ \hat{d}(\nu,f) &= \frac{1}{\Delta x \Delta t} \sum^\infty_{n=-\infty}\delta(\nu - \frac{n}{\Delta x}) \sum^\infty_{n=-\infty}\delta(f - \frac{n}{\Delta t}) \end{align} $$Thus, in total, we have the final analytical expression for the two-dimensional fourier transform of the signal:
$$ \begin{align} \hat{s}(\nu,f) &= \frac{A i}{8 \Delta x \Delta t} \\ & \times\Big[ \sum^\infty_{n=-\infty}\Big(L sinc(L(\nu - \frac{n}{\Delta x} + k)) + \frac{1}{2} L sinc(L(\nu - \frac{n}{\Delta x} + k) - 1) + \frac{1}{2} L sinc(L(\nu - \frac{n}{\Delta x} + k) + 1) \Big)\\ &\times \sum^\infty_{n=-\infty}\Big(T sinc(T(f - \frac{n}{\Delta t} + \omega)) + \frac{1}{2} T sinc(T(f - \frac{n}{\Delta t} + \omega) - 1) + \frac{1}{2} T sinc(T(f - \frac{n}{\Delta t} + \omega) + 1) \Big) \\ &- \sum^\infty_{n=-\infty}\Big(L sinc(L(\nu - \frac{n}{\Delta x} - k)) + \frac{1}{2} L sinc(L(\nu - \frac{n}{\Delta x} - k) - 1) + \frac{1}{2} L sinc(L(\nu - \frac{n}{\Delta x} - k) + 1) \Big)\\ &\times \sum^\infty_{n=-\infty}\Big(T sinc(T(f - \frac{n}{\Delta t} - \omega)) + \frac{1}{2} T sinc(T(f - \frac{n}{\Delta t} - \omega) - 1) + \frac{1}{2} T sinc(T(f - \frac{n}{\Delta t} - \omega) + 1) \Big) \Big] \\ \end{align} $$The function Jason2_py.model.power_spectrum() computes the above result over a range of $n = [-50, 50]$. Let's use this to model the principal semi-diurnal tides, M2 (lunar) and S2 (solar).
# M2 and S2 tidal periods (hours)
periods = [12.4206012,12]
semidiurnals = ['M2 model', 'S2 model']
for t_num, t_name in enumerate(tracks):
# unpack parameters from data
track = models[t_name]['data']
dx = track['dx']
dt = track['dt']
nu_vector = track['wavenumbers']
f_vector = track['frequencies']
L = dx * len(nu_vector)
T = dt * len(f_vector)
fft_2D = track['fft 2D']
for i,tide in enumerate(semidiurnals):
model = models[t_name][tide]
k = model['wavenumber']
# find alias frq from Jason-2 repeat period
w = Jason2.helper.find_alias_frq(9.9156,periods[i]/24)
model['alias frq'] = w
# simulate spectrum, unit amplitude
model_spectrum = Jason2.model.power_spectrum(1,L,nu_vector,k,T,w,f_vector,dx,dt)
# fit model to predicted tidal peak, save values for plotting
model['fft 2D'] = Jason2.model.fit_models([[k,w]],nu_vector,f_vector,fft_2D, [model_spectrum]) * model_spectrum
model['wavenumbers'] = nu_vector
model['frequencies'] = f_vector
model['frequencies'] = f_vector
np.seterr(divide='ignore')
# put plots into dict
_model_p = dict()
for i, track in enumerate(tracks):
for tide in ['M2 model', 'S2 model']:
_model_p[(tide.replace(" model",""),'Track ' + str(i + 1))] = pf.fft_2D_plt(models[track][tide],'model')
# make drop-down list
model_p = hv.HoloMap(_model_p, kdims=['Tide','Track'])
# plot overlay
(track_plot +
power_plot.opts(width =300,
height = 300,
title='data',
colorbar=False) +
model_p.opts(width =300,
height = 300,
title='model')).opts(title='Figure 2')
# initialize dict to hold plots
_f = dict()
_nu = dict()
for i, track in enumerate(tracks):
# get data spectrum
data = models[track]['data']
for tide in ['M2 model', 'S2 model']:
# get model spectrum
model = models[track][tide]
val_f = model['alias frq']
val_nu = model['wavenumber']
# make individual line plots
mf = pf.line_plt(model,'f',val_nu).relabel('model')
df = pf.line_plt(data,'f',val_nu).relabel('data')
mnu = pf.line_plt(model,'nu',val_f).relabel('model')
dnu = pf.line_plt(data,'nu',val_f).relabel('data')
# make overlay dicts
_f[(tide.replace(" model",""),'Track ' + str(i + 1))] = (df * mf).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
_nu[(tide.replace(" model",""),'Track ' + str(i + 1))] = (dnu * mnu).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
# # Make drop down list
f = hv.HoloMap(_f, kdims=['Tide','Track']).opts(infer_projection=True)
nu = hv.HoloMap(_nu, kdims=['Tide','Track']).opts(infer_projection=True)
# # plot layout
(f + nu).opts(title='Figure 3',shared_axes=False)
Anisotropic Wave with Frequency Doppler Shift
Comparing our M2 model to the data in frequency space, we find that the tidal peak is blue-shifted and broader than the model predicts. A possible explanation for this discrepancy is the effect of Doppler shifting by mean flows South of the Hawaiian Ridge. To examine this, we refine our model to include a time-dependent advection term using data from Copernicus. The Dopper shifted frequencies are calculated as:
$$\omega_t = \omega + {\overline{U}} \cdot \nu$$where $\overline{U}$ is the current velocity averaged over the water column, $\nu$ is the tidal wavenumber, and $\omega$ is the alias frequency.
First, We can plot the surface velocities near the tracks (Figure 4):
_U = {'Track ' + str(i + 1) : pf.U_plt(current_data[track],models[track]['data'],"") for i, track in enumerate(tracks)}
U = hv.HoloMap(_U, kdims='Track').opts(infer_projection=True)
U.opts(title='Figure 4: Surface Current Magnitudes')
Thus, from the average surface velocity, we would expect the maximum magntiude of Doppler-shifting to be:
for i,track in enumerate(tracks):
theta = np.arctan(np.diff(models[track]['data']['lat range']) / np.diff(models[track]['data']['lat range']))
u = np.mean(current_data[track]['u'])
v = np.mean(current_data[track]['v'])
mag = np.sqrt(u**2 + v**2) * np.cos(theta)
print(mag[0] * 24 * 3600 / 1000 * models[track]['M2 model']['wavenumber'], ' 1/day for track', str(i+1))
However, a better approximation would be to take the depth-averaged velocities, which can be found by solving for the vertical modes, as derived below.
The horizontal velocities of internal tides for a rigid-lid, flat-bottom ocean are described by the following Sturm-Liouville boundary value problem (1).:
$$\frac{d}{dz}\Big(\frac{1}{N^2}\frac{dF}{dz}\Big) + \frac{1}{c^2}F = 0, \\ \frac{dF}{dz} = 0 \ \ \ \ \text{at} \ \ \ \ z = 0 \ \ \ \ \text{ and} \ \ \ \ z= -H$$where $H$ is ocean depth and $N$ is the buoyancy frequency. The eigenvalue solutions (vertical modes), $F_n$, have associated eigenvalues of $-1/c_n^2$, and the tidal wavneumber is calculated using the dispersion relationship:
$$k_n = \frac{(\omega^2 - f^2)^{1/2}}{c_n}$$where $f$ is the coriolis parameter and $\omega$ is the tidal frequency (e.g. reported here). $F_n$ are solved using the dedalus framework with stratification from ECCO version 4. Estimates of current velocity are done by assuming horizontal velocities are given by
$$U = A_1 F_1 + A_2 F_2, \\ U = 0 \ \ \ \ \text{at} \ \ \ \ z = -H \ \ \ \ \text{and} \ \ \ \ U = U_1 \ \ \ \ \text{at} \ \ \ \ z = 0$$where $U_1$ is the current velocity reported by Copernicus projected on the Jason-2 track. In our model, Jason2_py.model.doppler_power_spectrum(), we simulate a power spectrum for each Jason-2 pass using the $\omega_t$ from that day and report the average of all simulations.
for t_num, t_name in enumerate(tracks):
data = models[t_name]['data']
eigenvector = models[t_name]['eigenvectors'][2,:]
c_data = current_data[t_name]
for tide in ['M2 model', 'S2 model']:
print('simulating track ' + ''.join([str(t_num + 1), ' ', tide]))
model = models[t_name][tide]
# evaluate doppler model
args = (data['lat range'],
data['lon range'],
data['frequencies'],
data['wavenumbers'],
model['alias frq'],
data['dx'] * len(data['wavenumbers']),
model['wavenumber'],
data['dt'] * len(data['frequencies']),
data['dx'],
data['dt'],
data['fft 2D'])
w_t,power_spec = Jason2.model.doppler_power_spectrum(eigenvector,c_data,Hawaii_SA[t_name],args)
# fit model to tidal wavenumber and alias frequency
coeff = Jason2.model.fit_models([[args[6],np.mean(w_t)]],args[3],args[2],args[-1], [power_spec])
fft_2D = coeff * power_spec
# store values
models[t_name]['doppler ' + tide] = {'wavenumber' : model['wavenumber'],
'alias frq' : model['alias frq'],
'shifted frequencies' : w_t,
'fft 2D' : fft_2D,
'wavenumbers' : data['wavenumbers'],
'frequencies' : data['frequencies']}
_wt = dict()
for i, track in enumerate(tracks):
frequencies, edges = np.histogram(models[track]['doppler M2 model']['shifted frequencies'], 20)
title = 'Alias Frequency: %.5f 1/day'%models[track]['M2 model']['alias frq']
_wt['Track ' + str(i + 1)] = hv.Histogram((edges, frequencies)).opts(xlabel = '𝜔ₜ (1/day)',
title = title,
fontsize = {'labels':12})
wt = hv.HoloMap(_wt, kdims='Track').opts(title='Figure 5')
wt
# put plots into dict
_doppler = dict()
for i, track in enumerate(tracks):
for tide in ['M2 model', 'S2 model']:
_doppler[(tide.replace(" model",""),'Track ' + str(i + 1))] = pf.fft_2D_plt(models[track]['doppler ' + tide],'doppler model')
# make drop-down list
doppler = hv.HoloMap(_doppler, kdims=['Tide','Track']).opts(width =300,height = 300)
# plot overlay
p1 = power_plot.opts(width = 250, height = 300,title='data', colorbar=False)
p2 = model_p.opts(width = 250,height = 300,title='model',colorbar=False)
(p1 + p2 + doppler).opts(title='Figure 6')
# initialize dict to hold plots
_frq = dict()
_om = dict()
for i, track in enumerate(tracks):
# get data spectrum
data = models[track]['data']
for tide in ['M2 model', 'S2 model']:
# get model spectrum
model = models[track][tide]
dmodel = models[track]['doppler ' + tide]
val_f = model['alias frq']
val_nu = model['wavenumber']
# make individual line plots
df = pf.line_plt(data,'f',val_nu).relabel('data')
mf = pf.line_plt(model,'f',val_nu).relabel('model')
dmf = pf.line_plt(dmodel,'f',val_nu).relabel('doppler model')
dnu = pf.line_plt(data,'nu',val_f).relabel('data')
mnu = pf.line_plt(model,'nu',val_f).relabel('model')
dmnu = pf.line_plt(dmodel,'nu',val_f).relabel('doppler model')
# make overlay dicts
_frq[(tide.replace(' model',''),'Track ' + str(i + 1))] = (df * mf * dmf).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
_om[(tide.replace(' model',''),'Track ' + str(i + 1))] = (dnu * mnu * dmnu).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
# Make drop down list
frq = hv.HoloMap(_frq, kdims=['Tide','Track']).opts(infer_projection=True)
om = hv.HoloMap(_om, kdims=['Tide','Track']).opts(infer_projection=True)
# plot layout
(frq.opts(width=450) + om.opts(width=400)).opts(title='Figure 7',shared_axes=False)
Below, we have a summary of mean square error:
# ranges (track x frequency range x wavenumber range)
ranges = np.array([[[.008,.03],[0.0055,0.008]],
[[.01,.03],[.003,.0085]],
[[0.000,.04,],[.005,.0075]],
[[.01,.025],[.005,.009]]])
# calculate mean square error
mse = np.zeros((4,2,2))
_range = dict()
for t_num,t_name in enumerate(tracks):
# get ranges
f_range = ranges[t_num,0,:].tolist()
nu_range = ranges[t_num,1,:].tolist()
_range['Track ' + str(t_num + 1)] = pf.fft_2D_plt(models[t_name]['data'],'').opts(xlim=tuple(nu_range),
ylim=tuple(f_range))
# get data spectrum
power_spectrum = models[t_name]['data']
f_vector = power_spectrum['frequencies']
nu_vector = power_spectrum['wavenumbers']
for tide_num, tide in enumerate(['M2 model', 'S2 model']):
# get model specta
model = models[t_name][tide]['fft 2D']
dmodel = models[t_name]['doppler ' + tide]['fft 2D']
# solve for mean square error
mse[t_num,tide_num,0] = Jason2.find_square_error(f_range,
nu_range,
f_vector,
nu_vector,
power_spectrum['fft 2D'],
model)
mse[t_num,tide_num,1]= Jason2.find_square_error(f_range,
nu_range,
f_vector,
nu_vector,
power_spectrum['fft 2D'],
dmodel)
# output
range_plt = hv.HoloMap(_range, kdims='Track').opts(framewise=True,
width=400,height=300)
d = dict()
d['M2'] = pd.DataFrame(columns=['Model', 'Doppler Model',],
data=mse[:,0,:],
index = ['Track %.f'%(i+1) for i in range(4)])
d['S2'] = pd.DataFrame(columns=['Model', 'Doppler Model',],
data=mse[:,1,:],
index = ['Track %.f'%(i+1) for i in range(4)])
pd.set_option('display.float_format', '{:.2E}'.format)
df = pd.concat(d, axis=1,names=['Tidal Constituent:'])
print('Table 1: Mean Square Error')
display(HTML(df.to_html()))
range_plt.opts(title='Figure 8: Semi-Diurnal Peaks')
Thus over the ranges selected, we see that the Doppler Model reduces the mean square error for all but track 4. Let's now see if we can pick up the same tidal signals from the descending tracks.
Load Jason-2 Data
# get existing cycles:
f = open("cycles.txt", "r")
cycles = [line.strip()[1:4] for line in f]
# define reference cycle
ref_cycle = '100'
# target start point and track length
starts = [[22.450964, -168.03407],
[21.74352, -164.891574],
[20.98662, -161.731304],
[18.619012, -160.738939]]
track_length = 1000
direction = 'S'
# download and read ssha data for ascending tracks South of Hawaiian Ridge
Hawaii_SD_raw = Jason2.nested_dict()
for track in range(4):
# define start and reference file:
start = starts[track]
reference_file = Jason2.download.find_track(start,ref_cycle,direction = direction)
track_name = 'SD' + str(track + 1)
# download, read data from each cycle
for cycle in cycles:
cycle_data = Jason2.download.read_track(start,track_length,ref_cycle,cycle,reference_file)
if cycle_data == None:
continue
Hawaii_SD_raw[track_name][cycle] = cycle_data
Load Copernicus Data (set directory name as necessary)
# set directories to retrieve nc files, initialize data structure
dir_names = ['D:/copernicus_current_data/{year}/{year}'.format(year=year) for year in np.arange(2008,2017)]
currents_SD = dict()
Hawaii_SD = Jason2.nested_dict()
lat_range = np.zeros((4,2))
lon_range = np.zeros((4,2))
for t_num, t_name in enumerate(Hawaii_SD_raw.keys()):
# preprocess data
Hawaii_SD[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SD_raw[t_name],starts[t_num],masked_tol=5)
# find latitude and longitude range
lats = np.array([[Hawaii_SD[t_name][cycle]['lat'][0],Hawaii_SD[t_name][cycle]['lat'][-1]] for cycle in Hawaii_SD[t_name].keys()])
lons = np.array([[Hawaii_SD[t_name][cycle]['lon'][0],Hawaii_SD[t_name][cycle]['lon'][-1]] for cycle in Hawaii_SD[t_name].keys()])
lat_range = [min(lats[:,1]), max(lats[:,0])]
lon_range = [min(lons[:,0]), max(lons[:,1])]
# download datasets from .nc directory (change directory name as necessary)
time,lat,lon,u,v = Jason2.download.read_current_data(lat_range,lon_range,dir_names)
currents_SD[t_name] = dict(zip(['time','lat','lon','u','v'],[time,lat,lon,u,v]))
# Jason-2
Hawaii_SD_raw = np.load('./test_data/Hawaii_SD_raw.npy',allow_pickle='TRUE').item()
# Copernicus
current_data = np.load('./test_data/current_dataSD.npy',allow_pickle='TRUE').item()
# ECCOv4
tracksd = ['SD' + str(i+1) for i in range(4)]
for track in tracksd:
modes = 1 / np.genfromtxt('./test_data\eigenvectors\modes%s.csv'%track,delimiter = ',')[:,1]
models[track]['M2 model']['wavenumber'] = modes[0]
models[track]['S2 model']['wavenumber'] = modes[1]
models[track]['eigenvectors'] = np.genfromtxt('./test_data\eigenvectors\%s.csv'%track,delimiter = ',')
Unlike the ascending tracks, we assume that waves do not propagate collinearly with the Jason-2 tracks. Thus, the effective wavenumbers used to describe the observed 2D power spectra are given as:
$$k_{eff} = k \ cos(\theta)$$where $k$ is the expected tidal wavenumber, and $\theta$ is the angle between the tide beam and the track. Since we wnat to model the strong semi-dirunal signals found in the ascending tracks, we have chosen 4 descending tracks whose midpoints are closest to the respective midpoints of the 4 ascending tracks. This choice was made because the Hann window places greater weight on values near the center of each track.
Below, we calculate $\theta$ for each track by simply taking the angle between each ascending-descending track pair:
Hawaii_SD = Jason2.nested_dict()
_trkd = dict()
_trk = dict()
tracksd = ['SD' + str(i+1) for i in range(4)]
starts = [[22.450964, -168.03407],
[21.74352, -164.891574],
[20.98662, -161.731304],
[18.619012, -160.738939]]
thetas = [0] * 4
for t_num, t_name in enumerate(tracksd):
# preprocess data
Hawaii_SD[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SD_raw[t_name],starts[t_num],masked_tol=5)
# find average lat and lon range
lats = np.array([[Hawaii_SD[t_name][cycle]['lat'][0],Hawaii_SD[t_name][cycle]['lat'][-1]] for cycle in Hawaii_SD[t_name].keys()])
lons = np.array([[Hawaii_SD[t_name][cycle]['lon'][0],Hawaii_SD[t_name][cycle]['lon'][-1]] for cycle in Hawaii_SD[t_name].keys()])
lat_range = [np.mean(lats[:,0]), np.mean(lats[:,1])]
lon_range = [np.mean(lons[:,0]), np.mean(lons[:,1])]
# save info to models
models[t_name]['data'].update({'lat range' : lat_range,
'lon range' : lon_range})
# find track angles:
data_a = models['SA' + str(t_num + 1)]['data']
data_d = models[t_name]['data']
v_a = [np.diff(data_a['lat range'])[0],np.diff(data_a['lon range'])[0]]
v_d = [np.diff(data_d['lat range'])[0],np.diff(data_d['lon range'])[0]]
thetas[t_num] = np.arccos(np.dot(v_d/np.linalg.norm(v_d),v_a/np.linalg.norm(v_a)))
# create plot dicts
_trkd['Track ' + str(t_num+1)] = pf.track_plt(models[t_name]['data'],"Track %.f"%(t_num+1))
_trk['Track ' + str(t_num+1)] = (track_dict['Track ' + str(t_num+1)] * _trkd['Track ' + str(t_num+1)]).opts(projection = ccrs.PlateCarree(),
title='Figure 9: θ = %.3f'%thetas[t_num])
# make drop down plots
trk = hv.HoloMap(_trk,kdims=['Track'])
trk
Plotting the 2D power spectra of the descending tracks:
_powerd = dict()
for t_num, t_name in enumerate(tracksd):
# preprocess data
Hawaii_SD[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SD_raw[t_name],starts[t_num],masked_tol=5)
# find power spectrum
fft_2D, dx, dt, nu_vector, f_vector = Jason2.find_2D_power_spectrum(Hawaii_SD[t_name])
# save info to models
models[t_name]['data'].update({'fft 2D' : fft_2D,
'dx': dx,
'dt': dt,
'wavenumbers' : nu_vector,
'frequencies' : f_vector})
# Make heatmaps
_powerd['Track ' + str(t_num + 1)] = pf.fft_2D_plt(models[t_name]['data'], '')
# make drop-down plots
powerd = hv.HoloMap(_powerd, kdims=['Track'])
(trk.opts(title='',infer_projection=True,height=200,width=200) +
power_plot.opts(title='Ascending',height=300,width=300) +
powerd.opts(title='Descending',height=300,width=300)).opts(title='Figure 10: 2D Power Spectra')
As expected, we are fitting to the peaks where the tidal wavenumber and alias frequency are opposite in sign. Now, let's turn to modeling. Plotting the surface currents we see similar magnitudes as for the ascending tracks, so we expect the Doppler shifting of the alias frequency to be a factor.
_Ud = {'Track ' + str(i + 1) : pf.U_plt(current_data[track],models[track]['data'],"") for i, track in enumerate(tracksd)}
Ud = hv.HoloMap(_Ud, kdims='Track').opts(infer_projection=True)
Ud.opts(title='Figure 11: Mean Flow at Surface')
Let's compute both the simple and Doppler model, and update our dictionaries:
# M2 and S2 tidal periods (hours)
periods = [12.4206012,12]
semidiurnals = ['M2 model', 'S2 model']
for t_num, t_name in enumerate(tracksd):
# get data
c_data = current_data[t_name]
data = models[t_name]['data']
eigenvector = models[t_name]['eigenvectors'][2,:] # 2nd vertical mode
dx = data['dx']
dt = data['dt']
nu_vector = data['wavenumbers']
f_vector = data['frequencies']
L = dx * len(nu_vector)
T = dt * len(f_vector)
fft_2D = data['fft 2D']
for i,sdiurnal in enumerate(semidiurnals):
model = models[t_name][sdiurnal]
k = model['wavenumber']
# multiply k by track angle
k_eff = k * np.cos(thetas[t_num])
# find alias frq from Jason-2 repeat period
w = Jason2.helper.find_alias_frq(9.9156,periods[i]/24)
model['alias frq'] = w
#-------Model 1: No Doppler shift--------
print('simulating track ' + ''.join([str(t_num + 1), ' ', sdiurnal]))
model_spectrum = Jason2.model.power_spectrum(1,L,nu_vector,k_eff,T,w,f_vector,dx,dt)
model['fft 2D'] = Jason2.model.fit_models([[k_eff,w]],nu_vector,f_vector,fft_2D, [model_spectrum]) * model_spectrum
model['wavenumbers'] = nu_vector
model['frequencies'] = f_vector
model['frequencies'] = f_vector
model['effective wavenumber'] = k_eff
#-------Model 2: Doppler Shift---------
args = (data['lat range'],
data['lon range'],
data['frequencies'],
data['wavenumbers'],
model['alias frq'],
data['dx'] * len(data['wavenumbers']),
k_eff,
data['dt'] * len(data['frequencies']),
data['dx'],
data['dt'],
data['fft 2D'])
w_t,power_spec = Jason2.model.doppler_power_spectrum(eigenvector,c_data,Hawaii_SD[t_name],args)
# fit model to tidal wavenumber and average Doppler-shifted frequency
coeff = Jason2.model.fit_models([[args[6],np.mean(w_t)]],args[3],args[2],args[-1], [power_spec])
fft_2D = coeff * power_spec
models[t_name]['doppler ' + sdiurnal] = {'wavenumber' : model['wavenumber'],
'effective wavenumber' : k_eff,
'alias frq' : model['alias frq'],
'shifted frequencies' : w_t,
'fft 2D' : fft_2D,
'wavenumbers' : data['wavenumbers'],
'frequencies' : data['frequencies']}
# put plots into dict
_modeld = dict()
_dmodeld = dict()
for i, track in enumerate(tracksd):
for tide in ['M2 model', 'S2 model']:
_modeld[(tide.replace(" model",""),'Track ' + str(i + 1))] = pf.fft_2D_plt(models[track][tide],'model')
_dmodeld[(tide.replace(" model",""),'Track ' + str(i + 1))] = pf.fft_2D_plt(models[track]['doppler ' + tide],'model')
# make drop-down list
modeld = hv.HoloMap(_modeld, kdims=['Tide','Track'])
dmodeld = hv.HoloMap(_dmodeld, kdims=['Tide','Track'])
trkd = hv.HoloMap(_trkd, kdims =['Track']).opts(infer_projection=True)
# plot overlay
(powerd.opts(width =250,
height = 300,
title='data',
colorbar=False) +
modeld.opts(width =250,
height = 300,
title='model',
colorbar=False) +
dmodeld.opts(width =300,
height = 300,
title='Doppler model')).opts(title='Figure 12: 2D models')
(powerd + modeld).opts(shared_axes=True)
# initialize dict to hold plots
_frqd = dict()
_omd = dict()
for i, track in enumerate(tracksd):
# get data spectrum
data = models[track]['data']
for tide in ['M2 model', 'S2 model']:
# get model spectrum
model = models[track][tide]
dmodel = models[track]['doppler ' + tide]
val_f = model['alias frq']
val_nu = model['effective wavenumber']
# make individual line plots
dfd = pf.line_plt(data,'f',val_nu).relabel('data')
mfd = pf.line_plt(model,'f',val_nu).relabel('model')
dmfd = pf.line_plt(dmodel,'f',val_nu).relabel('doppler model')
dnud = pf.line_plt(data,'nu',val_f).relabel('data')
mnud = pf.line_plt(model,'nu',val_f).relabel('model')
dmnud = pf.line_plt(dmodel,'nu',val_f).relabel('doppler model')
# make overlay dicts
_frqd[(tide.replace(' model',''),'Track ' + str(i + 1))] = (dfd * mfd * dmfd).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
_omd[(tide.replace(' model',''),'Track ' + str(i + 1))] = (dnud * mnud * dmnud).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
# Make drop down list
frqd = hv.HoloMap(_frqd, kdims=['Tide','Track']).opts(infer_projection=True)
omd = hv.HoloMap(_omd, kdims=['Tide','Track']).opts(infer_projection=True)
# plot layout
(frqd.opts(width=450) + omd.opts(width=400)).opts(title='Figure 13',shared_axes=False)
Descending Track Models: Mean Square Error
# ranges (track x frequency range x wavenumber range)
rangesd = np.array([[[.01,.03],[-.008,-.0035]],
[[.01,.03],[-.008,-.0035]],
[[0.01,.03,],[-.006,-.003]],
[[.015,.025],[-.008,-.004]]])
# calculate mean square error
msed = np.zeros((4,2,2))
_ranged = dict()
for t_num,t_name in enumerate(tracksd):
# get ranges
f_range = rangesd[t_num,0,:].tolist()
nu_range = rangesd[t_num,1,:].tolist()
_ranged['Track ' + str(t_num + 1)] = pf.fft_2D_plt(models[t_name]['data'],'').opts(xlim=tuple(nu_range),
ylim=tuple(f_range))
# get data spectrum
power_spectrum = models[t_name]['data']
f_vector = power_spectrum['frequencies']
nu_vector = power_spectrum['wavenumbers']
for tide_num, tide in enumerate(['M2 model', 'S2 model']):
# get model specta
model = models[t_name][tide]['fft 2D']
dmodel = models[t_name]['doppler ' + tide]['fft 2D']
# solve for mean square error
msed[t_num,tide_num,0] = Jason2.find_square_error(f_range,
nu_range,
f_vector,
nu_vector,
power_spectrum['fft 2D'],
model)
msed[t_num,tide_num,1]= Jason2.find_square_error(f_range,
nu_range,
f_vector,
nu_vector,
power_spectrum['fft 2D'],
dmodel)
# output
range_pltd = hv.HoloMap(_ranged, kdims='Track').opts(framewise=True,
width=400,height=300)
dd = dict()
dd['M2'] = pd.DataFrame(columns=['Model', 'Doppler Model',],
data=msed[:,0,:],
index = ['Track %.f'%(i+1) for i in range(4)])
dd['S2'] = pd.DataFrame(columns=['Model', 'Doppler Model',],
data=msed[:,1,:],
index = ['Track %.f'%(i+1) for i in range(4)])
pd.set_option('display.float_format', '{:.2E}'.format)
dfd = pd.concat(dd, axis=1,names=['Tidal Constituent:'])
print('Table 2: Mean Square Error')
display(HTML(dfd.to_html()))
range_pltd.opts(title='Figure 8: Semi-Diurnal Peaks')